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ABSTRACT 


The  SIGMA  package  is  a  set  of  FORTRAN  subprograms,  using  double¬ 
precision  floating-point  arithmetics,  which  attempts  to  find  a  global 
minimizer  of  a  real-valued  function  f(x)  =  f(x-j,  ...  ,  x^)  of  N 
real  variables  x^ ,  ...  ,  x^. 


AMS  (MOS)  Subject  Classifications:  65K10,  60H10,  49D25 

Key  Words:  Algorithms,  Global  Optimization,  Stochastic  Differential 
Equations  , 

Work  Unit  Number  5  -  Optimization  and  Large  Scale  Systems 


Dipartimento  di  Matematica,  Universita  di  Bari,  70125  Bari  (Italy). 

**Istituto  di  Fisica,  2a  Universita  di  Roma  "Tor  Vergata",  Via  Orazio 
Rairaondo,  00173  (La  Romanina)  Roma  (Italy). 

***Istituto  di  Matematica,  Universita  di  Salerno,  84100  Salerno  (Italy). 


Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-80-C-0041 . 

The  research  reported  in  this  paper  has  been  made  possible  through  the 
support  and  sponsorship  of  the  U.S.  Government  through  its  European  Research 
Office  of  the  U.S.  Army  under  Contract  No.  DAJA-37-81-C-0740  with  the 
University  of  Camerino,  Italy. 

/I 


SIGNIFICANCE  AND  EXPLANATION 

The  paper  gives  a  detailed  description  of  a  FORTRAN  IV  program  based  on 
a  new  method  for  finding  a  "global"  (or  "absolute")  minimizer  of  a  function 
of  N  real  variables,  i.e.  the  point  x  in  N-dimensional  space  (or  possibly 
one  of  the  points)  such  that  not  only  the  function  Increases  if  one  moves 
away  from  x  in  any  direction,  ("local"  or  "relative"  minimum),  but  also  such 
that  no  other  point  exists  where  f  has  a  lower  value. 

The  method,  which  was  first  proposed  by  the  present  authors  in  a  paper 
which  is  to  appear  in  the  Journal  of  Optimization  Theory  and  Applications,  is 
based  on  ideas  from  statistical  mechanics,  and  looks  for  a  point  of  global 
minimum  by  following  the  solution  trajectories  of  a  stochastic  differential 
equation  representing  the  motion  of  particle  (in  N-space)  under  the  action 
of  a  potential  field  and  of  a  random  perturbing  force.  ^ 

The  complete  algorithm  is  described  in  detail iff  the  companion  paper  by 
the  same  authors:  A  Global  Optimization  Algopftfun  using  Stochastic  Differ¬ 
ential  Equations,  which  also  contai_ns~*"summary  of  the  results  of  extensive 
numerical  testing. 

^The  tests  were  performed  by  running  the  program  on  an  extensive  set  of 
carefully  selected  test  problems  of  varying  difficulty,  and  the  performance 
was  remarkably  successful,  even  on  very  hard  problems  (e.g.  problems  with  a 
single  point  of  global  minimum  and  up  to  about  lO^  points  of  non-global 
minimum).  /o  -y®  ib  ic tk  fwea 

It  is  finally  to  be  noted  that  the  majority  of  the  optimization  methods 
presently  available  deal  with  the  local  optimization  problem,  and  that  no 
methods  of  comparable  power  seem  to  be  available  in  the  field  of  global 
optimization. 

The  program  is  essentially  designed  for  portability  on  different  types 
of  computers,  since  it  meets  the  specifications  of  PFORT,  a  portable  subset 
of  ANS  FORTRAN  developed  at  Bell  Laboratories. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  HRC,  and  not  with  the  authors  of  this  report. 
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1.  METHOD 

The  algorithm  used  by  SIGM  is  described  in  detail  in  ref.  [1] . 

A  global  minimizer  of  f(x)  is  sought  by  monitoring  the  values 
of  f(x)  along  trajectories  generated  by  a  suitable  discretization 
of  the  stochastic  differential  equation 
d£  *  -Vf(Odt  +  e(t)dw 
with  initial  condition: 

£(0)  -  Xo 

where  Vf  is  the  gradient  of  f,  w(t)  is  an  N-dimensional  standard 
Wiener  process,  and  the  "noise  coefficient"  e(t)  is  a  positive  func¬ 
tion.  The  discretization  has  the  form 

"Sk  '  *  c<V-*k  V  k* 0-1*2"  ••• 

So  •  & 

where  h^  is  the  time  integration  step  length,  S  7(5^)  is  computed  as 
a  finite-differences  approximation  to  the  directional  derivative  of  f 
in  a  randomly  chosen  direction,  and  is  a  random  sample  from  an  N- 
dimensional  standard  gaussian  distribution. 
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We  consider  the  simultaneous  evolution  of  a  number  of  tra¬ 

jectories  during  an  "observation  period"  having  the  duration  of  a  given 
number  of  time  integration  steps,  and  within  which  the  noise  coeffi¬ 

cient  e(t)  of  each  trajectory  is  kept  at  a  constant  value  £p,  while 

m 

are  automatically  adjusted  for  each  trajectory  by  the  algorithm. 

At  the  end  of  every  observation  period  a  comparison  is  made  between 
the  trajectories:  one  of  the  trajectories  is  discarded,  all  other  tra¬ 
jectories  are  naturally  continued  in  the  next  observation  period,  and  one 
of  them  is  selected  for  "branching",  that  is  for  generating  also  a  second 
continuation  trajectory  which  differs  from  the  first  one  only  in  the  start¬ 
ing  values  for  and  Ax^,  and  is  considered  as  having  the  same  "past 
history"  of  the  first. 

The  nunber  of  N^j  of  simultaneously  evolving  trajectories  remains 
therefore  unaffected,  and  the  second  continuation  trajectory  takes  the 
place,  from  a  program- implementation  point  of  view,  of  the  discarded  tra¬ 
jectory. 

The  set  of  simultaneous  trajectories  is  considered  as  a  single 
trial,  and  the  complete  algorithm  is  a  set  of  repeated  trials.  A  single 
trial  is  stopped,  at  the  end  of  an  observation  period,  if  a  maximum  given 
nunber  of  observation  periods  has  been  reached,  or  if  all  the 

final  values  of  f(x)  (except  for  the  discarded  trajectory)  are  equal 
(\%ithin  numerical  tolerances,  and  possibly  at  different  points  x)  to 
their  minimum  value  ("unifonn  stop"  at  the  level  .  In 

the  former  case  the  trial  is  considered  unsuccessful,  while  in  the  latter 
case  a  comparison  is  made  between  the  common  final  function  value 


the  steplength  h^  and  the  spatial  increment  Ax^  for  confuting 
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and  the  current  best  minimum  function  value  fgp-p  found  so  far  from 
algorithm  start:  if  >  f0pr  the  trial  is  a8ain  considered  unsucces 

ful;  and  if  *  ^OPT  (within  nunerical  tolerances)  the  trial  is 

considered  successful  at  the  level  fQpp. 

The  trials  are  repeated  with  different  operating  conditions  (initial 
point  Xq,  maximum  trial  length  »  see<*  t^ie  n°ise  generator, 

policy  for  selecting  the  starting  value  for  ep  in  the  second  continua¬ 
tion  trajectory  after  branching,  and  trial-start  value  for  ep)  and  the 
complete  algorithm  is  stopped  —  at  the  end  of  a  trial  —  if  a  given 
number  NsuC  of  uniform  stops  at  the  current  fQpj.  level  has  been  ob¬ 
tained,  or  if  a  given  maximum  number  of  trials  has  been  reached: 

success  of  the  algorithm  is  claimed  if  at  least  one  uniform  stop  occurred 
at  the  final  value  of  f^pj. 


2.  DESCRIPTION  OF  THE  PACKAGE 


The  algorithm  used  by  SIQ4A  (see  sect.  2  and  ref.  [1])  has  been 
coded  in  the  form  of  a  set  of  FORTRAN  subprograms,  using  double -precis ion 
floating-point  arithmetics,  which  are  described  below. 

2.1.  Language 

All  the  coding  is  written  in  FORTRAN  IV  and  meets  the  specifications 
of  PFORT,  a  portable  subset  of  A.N.S.  FORTRAN  (ref.  [2]).  The  FORTRAN 
implicit  type  definition  for  integers  is  used  throughout;  all  non-integer 
variables  are  double  precision. 

2.2.  Description  of  the  Subprograms 

The  SICMA  package  consists  of  a  principal  subroutine  SIGMA,  a  set 
of  27  auxiliary  subroutines,  INIT,  REINIT,  TRIAL,  GENEVA,  PERIOD,  BRASI, 
ORDER,  CCMPAS,  STEP,  SSTEP,  NEWH,  DERFOR,  DERCEN,  RCLOPT,  STOOPT,  RANGE, 
INISCA ,  NOSCA,  SEGSCA,  VARSCA,  CIWSCA,  ACTSCA,  MOVSCA,  UPDSCA,  ALKMJT, 
GAUSRV,  UNITRV;  a  set  of  7  auxiliary  functions,  IPREC,  IPRECE,  FUNCT0, 
ITOLCH,  EIGSCA,  CHAOS,  UNIFRN;  and  a  driver  subroutine  SI  (MAI  calling 
SIGMA.  The  subprograms  are  described  below,  The  user  interested  only  in 
the  use  of  SIGMA  may  jump  to  Section  3. 

We  may  group  the  subprograms  as  follows. 

a)  Subprograms  for  the  numerical  integration:  STEP,  SSTEP,  DERFOR, 
DERCEN,  FUNCT0,  RANGE,  NEWH.  The  value  of  the  function  f(x)  is  computed 
—  whenever  required  in  the  numerical  integration  process  —  by  calling  the 
function  FINCT0.  FUNCT0  rescales  the  variables  by  calling  VARSCA  (see  d)), 
calls  RANGE  to  take  care  of  the  cases  where  the  current  point  x  is 
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outside  the  admissible  range  ([1],  sect.  3.2.11)  calls  the  user-supplied 
function  FUNCT  (sect.  3.5.1)  to  compute  f  (x) ,  and  possibly  updates  the 
best  current  function  minimum  fg^  and  the  corresponding  minimizer  Xqpt 
by  calling  STOOPT  (see  c)).  The  basic  step  of  the  numerical  integration 
is  performed  by  SSTEP  which  calls  FUNCT0  to  compute  the  value  of  f (x) , 
and  UN'ITRV  (see  e))  to  compute  the  random  direction  along  which  the  direc¬ 
tional  derivative  is  to  be  computed  (see  [1],  sect.  3.2.2);  the  directional 
derivatives  are  computed  numerically  by  SSTEP,  with  forward  or  central 
finite  differences,  by  calling  DERFOR  or  DERCEN,  which  call  FUNCT0;  the 
first  half-step  ([1] ,  sect.  3.2.1)  is  accepted  or  rejected  ([1],  sect.  3.2.3) 
by  calling  NEWH  which  also  provides  the  updated  value  of  the  time  inte¬ 
gration  steplength  h^;  SSTEP  also  updates  the  cumulated  scaling  data  ([1] , 
sect.  3.2.12)  by  calling  CUMSCA  (see  d)),  and  updates  the  spatial  discreti¬ 
zation  increment  Ax^  based  on  the  results  of  calling  ITOLCH.  The  second 
half-step  ([1],  sect.  3.2.1)  is  performed  by  SSTEP  by  calling  GAUSRV 
(see  e))  and  can  be  accepted  or  rejected  ([1],  sect.  3.2.3).  The  subrou¬ 
tine  STEP  performs  the  single  integration  step  for  each  one  of  the  simul¬ 
taneous  trajectories  by  repeatedly  calling  SSTEP. 

b)  Subprograms  for  the  selection  of  the  trajectories:  BRASI,  ORDER, 
IPREC,  IPRECE,  CCMPAS.  The  selection  process  for  the  trajectories  ([1]  , 
sect.  3.2.6)  is  performed  by  the  subroutine  BRASI.  BRASI  first  updates 
the  trajectory  data  corresponding  to  the  elapsed  observation  period,  and 
then  asks  for  an  ordering  of  the  trajectories  by  calling  ORDER.  ORDER 
obtains  the  ordering  by  comparing  two  trajectories  on  the  basis  of  the 
past  history,  (by  calling  IPREC),  and  of  the  value  of  the  noise  coefficient 
Cp  (by  calling  IPRECE)  ([1],  sect.  3.2.6).  Based  on  the  ordering  provided 
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by  ORDER,  BRASI 

1)  discards  one  of  the  trajectories 

2)  performs  a  branching  on  another  trajectory,  i.e.  the 
trajectory  to  be  branched  gives  rise  to  two  "continuation" 
trajectories:  the  first  one  is  unperturbed,  and  the  second 

one  has  modified  values  for  e  and  for  the  initial  Ax,  ; 

P  * 

the  modified  values  are  obtained  from  the  old  ones  by  means 
of  random  multiplicative  factors  which  are  computed  with  the 
aid  of  random  number  generator  function  CHAOS  (see  e)). 

Since,  from  a  program  implementation  point  of  view*,  the  new  trajectory 
is  "moved"  in  the  "position"  of  the  discarded  one,  all  the  trajectory 
parameters  must  be  moved  to  the  new  position.  This  is  performed  directly 
by  BRASI  for  all  the  trajectory  data,  except  for  the  scaling  data  which 
are  moved  by  MOVSCA  (see  d)).  Finally  BRASI  calls  CQMPAS  in  order  to  examine 
the  stored  data  about  past  trajectories  from  the  point  of  view  of  their 
utility  to  the  only  user  of  such  data,  which  is  the  subroutine  IPREC,  and 
irrelevant  data  are  discarded. 

c)  Subprograms  for  general  management  of  the  complete  algorithm: 

SICMA,  INIT,  REINIT,  TRIAL,  GENEVA,  PERIOD,  ITOLCH,  RCLOPT,  STOOPT. 

GENEVA  performs  the  generation  of  the  set  of  trajectory  segments 
corresponding  to  the  current  observation  period  and  the  final  processing 
and  evaluation  of  the  trajectories.  GENEVA  first  updates  the  scaling 
arrays  containing  A  and  b  ([1],  sect.  3.2.12)  by  calling  SEGSCA  and 
UPDSCA  (see  d)) .  The  generation  of  the  trajectory  segments  is  performed 
by  GENEVA  by  calling  PERIOD. 


PERIOD  first  computes  the  duration  (in  accepted  steps)  of  the  observation 
period,  computes  all  the  integration  steps  by  repeatedly  calling  STEP  (see 
a))  and  finally  performs  the  trajectory  selection  by  calling  BR4SI  (see  b)). 

Finally  GENEVA  determines  some  end -of -segment  results  (FPFMIN,  FPFMAX, 
XPFMIN,  see  sect.  3.5.2)  using  the  rescaling  capabilities  of  SEGSCA  and 
VARSCA  (see  d)). 

The  subroutine  TRIAL  generates  a  trial  by  repeatedly  performing,  for 
every  observation  period, 

a  call  to  GENEVA  which  generates  the  simultaneous  trajectory 
segments,  and  performs  the  trajectory  selection, 
a  (possible)  call  to  PTSEG  which  performs  end-of- segment  output, 
a  check  of  the  (trial)  stopping  criteria,  with  the  aid  of  the 
function  ITOLCH, 

a  decision  about  activating  or  deactivating  the  scaling  of  the 
variables  (actions  performed  by  calling  ACTSCA  or  NOSCA) . 

The  subroutine  SICMA  is  the  principal  subroutine  of  the  package 
and  is  the  only  one  which  must  be  called  by  the  user  (apart  from  the  driver 
SIGMA1,  sect.  3 .4) . 

SICMA  manages  the  execution  of  the  complete  algorithm,  i.e.  of  a 
sequence  of  repeated  trials  performed  by  varying  a  number  of  operating 
conditions.  SIGMA  initializes  the  first  trial  by  calling  INIT,  and  the 
other  trials  by  calling  REINIT. 

For  each  trial  the  subroutine  SICMA 

enables  or  prevents  a  future  activation  (within  the  current  trial) 
of  the  scaling  of  the  variables  by  calling  INISCA  or  NOSCA 
actually  executes  the  trial  by  calling  TRIAL 
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updates  a  number  of  parameters  (using  ITOLCH  and  RCLOPT) 

checks  the  algorithm-stopping  criteria 

possibly  performs  end-of-trial  outputs  by  calling  PTSEG,  PTRIAL, 

and  PTKSUC  (see  3.5.2) 

The  subroutine  STOOPT  and  RCLOPT  respectively  "store"  and  "recall" 
the  current  values  of  the  best  minimum  FOPT  and  of  the  corresponding 
minimi zer  XOPT. 

d)  Subprograms  for  rescaling  the  variables:  INISCA,  NOSCA,  SEGSCA, 
VARSCA,  CUMSCA,  ACTSCA,  MDVSCA,  UPDSCA,  EIGSCA  ([1],  sect.  3.2.12). 

INISCA  initializes  the  coimion  area  /SCALE/  for  the  scaling  data. 

NOSCA  deactivates  the  rescaling. 

SEGSCA  selects  the  trajectory  which  must  be  rescaled. 

VARSCA  computes  the  rescaled  variables  Ax  +  b. 

CUMSCA  stores  cumulated  statistical  data  on  the  ill -conditioning 
of  f  (Ax  +  b ) . 

ACTSCA  activates  the  rescaling.  • 

MDVSCA  moves  the  scaling  data  from  the  first  to  the  second  con¬ 
tinuation  of  a  branched  trajectory. 

UPDSCA  updates  the  scaling  matrix  A  and  vector  b  by  calling 
EIGSCA  and  VARSCA. 

EIGSCA  computes  the  largest  eigenvalue  of  a  matrix  used  for  rescal¬ 
ing,  starting  from  randomly  chosen  estimates  (obtained  by  calling  UNITRV) 
of  the  corresponding  eigenvector. 

e)  Subprograms  for  pseudo-random  number  generation:  CHAOS,  UNIFRN, 
ALKNUT,  GAUSRV,  UNITRV. 
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CHAOS  generates  an  element  of  a  sequence  of  independent  pseudo  randan 
numbers,  each  one  having  one  out  of  four  possible  probability  distributions 
(standard  gaussian,  standard  Cauchy  ([lj,  sect.  3.2.7),  uniform  in  (-1,1) 
or  (0,1))  by  calling  UNIFRN. 

UNIFRN  generates  an  element  of  a  sequence  of  independent  pseudo¬ 
random  numbers  uniformly  distributed  in  (0,1),  by  calling  ALKNITT  and 
performing  a  further  (nonlinear)  randomization. 

ALKNUT  generates  an  element  of  a  sequence  of  independent  pseudo¬ 
random  numbers  (algorithm  of  Mitchell  and  Moore,  modified  as  suggested 
by  Brent,  see  ref.  [3]). 

GAUSRV  generates  an  element  of  a  sequence  of  independent  pseudo¬ 
random  N-vectors,  having  an  N-dimensional  standard  gaussian  probability 
distribution,  by  means  of  a  rejection  method,  and  based  on  uniformly 
distributed  (-1,1)  pseudo-random  numbers  obtained  by  calling  CHAOS. 

UNITRV  generates  an  element  of  a  sequence  of  independent  pseudo¬ 
random  N-vectors  uniformly  distributed  on  the  unit  sphere  in  RN. 


-9- 


1.  c 

2.  C  MAIN  PROGRAM  (SAMPLE  VERSION) 

3.  C  CALLS  SIGMA  VIA  THE  DRIVER  SUBROUTINE  SIGMA1 

4.  C 

5.  DOUBLE  PRECISION  X0,  XMIN,  FMIN 

6.  C 

7.  DIMENSION  X0(2),  XMIN(2) 

8  C 

9!  C  TEST  PROBLEM  DATA 

10.  C 

11.  C  PROBLEM  DIMENSION 

12.  N  =  2 

13.  C 

14!  C  INITIAL  POINT 

15.  X0(1)  -  M0 

16.  X0(2)  -  0.D0 

17.  C 

18.  C  SET  INPUT  PARAMETERS 

19.  NSUC  =  3 

20.  IPRINT  =  0 

21 .  C 

22*.  C  CALL  DRIVER  SUBROUTINE  SIGMA1 

23.  CALL  SIGMA1  (N,X0,NSUC,INPRINT,XMIN,FMIN,NFEV,IOUT) 

24.  C 

25.  STOP 

26.  END 


27.  DOUBLE  PRECISION  FUNCTION  FUNCT  (N,X) 

28  C 

29!  C  COMPUTES  THE  VALUE  AT  X  OF  THE  SIX-HUMP  CAMEL  FUNCTION 

30.  C 

31.  DOUBLE  PRECISION  X,XX,YY 

32.  DIMENSION  X(N) 

33.  XX  =  X(1  )*X(1 ) 

34.  YY  =  X ( 2 )*X ( 2 ) 

35.  FUNCT  =  ( (XX/3.D0-2. 1D0)*XX+4.D0)*XX+X(1 )*X(2) 

36.  1  +  4.D0*(YY-1.D0)*YY 

37.  RETURN 

38.  END 


Figure  1  --  List  of  the  Sample  Program 
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3.6.  Storage  Requirements 

The  SIGMA  package  contains  a  total  of  about  1900  statements  (includ¬ 
ing  some  700  coiment  lines) .  This  amounts  on  the  ASCII  FORTRAN  compiler 
(with  optimization  option)  of  the  UNIVAC  EXEC  8  operating  system  to  a  stor¬ 
age  requirement  of  about  4000  (36-bit)  words  for  the  instructions,  about 
3500  words  for  the  data,  and  about  14,000  words  for  the  COM40X  area.  The 
requirement  for  the  array  dimensions  are  4N  36 -bit  words. 

3.7.  Example 

T 

Let  N  «  2,  x  =  (x^^)  i  and  consider  the  six-humps  camel  func- 

1  6  4  2  4  2 

tion  f(x)  -  jX^  -  2.1  x^  +  4x.^  +  x^  +  4X2  -  4x2  which  has  four  non- 

T 

global  minima,  and  two  global  minima  at  x  ~  ±  (-0.089842,  0.71266)  where 
f  r  -1.03163.  The  sample  program  listed  in  fig.  1,  which  uses  the  easy- 
to-use  driver  SICMAl,  was  run  on  a  UNIVAC  1100/82  conputer  with  EXEC8 
operating  system  (level  38R5)  and  ASCII  FORTRAN  compiler  (version  10R1A) , 
starting  from  Xq  *  (0,0)^  and  with  NSUC  =  3. 

The  program  claimed  success  (IOUT  =  1)  stopping  correctly  at  one  of 
the  global  minimizers,  using  19660  function  evaluations.  The  printout  showed 
that  if  Ngyc  had  been  equal  to  1  (resp.  2) ,  the  minimum  would  have  been 
found  with  only  2697  (resp.  8343)  function  evaluations. 

ACKNOhlEDGEMENT :  One  of  us  (F.Z.)  gratefully  acknowledges  the  hospitality 
and  the  support  of  the  Mathematics  Research  Center  of  the  University  of 
Wisconsin  and  of  the  Department  of  Mathematical  Sciences  of  Rice  University 
where  part  of  this  work  was  done. 
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|ISTOP|  -  1  relative  difference  criterion  satisfied 
*  2  absolute  difference  criterion  satisfied 


•  3  both  criteria  satisfied 

The  sign  of  ISTOP  indicates  the  relationship  between  the  end-of-trial 
value  FTFMIN  and  the  best  current  minimum  value  FOPT  (which  is  updated 
whenever  a  function  value  is  computed) . 

ISTOP  >  0  FTFMIN  is  numerically  equal  (with  respect  to  at 

least  one  of  the  above  difference  criteria)  to  FOPT. 
ISTOP  <  0  FTFMIN  is  not  even  numerically  equal  to  FOPT  (and 

therefore  cannot  be  considered  an  acceptable  estimated 
global  minimum). 

ISTOPT  is  the  value  of  the  trial  stopping  indicator  ISTOP  correspond¬ 
ing  to  the  (current  or  past)  trial  where  FTFOPT  was  obtained,  with 
the  sign  which  is  updated  according  to  the  comparison  between  FTFOPT 
and  the  present  value  of  FOPT,  as  described  above.  The  final  value 
of  ISTOPT  is  returned  by  SIGMA,  as  the  value  of  the  output  indicator 
IOUT  (whenever  the  algorithm  was  started,  IOUT  +  -99,  see  above). 

The  subroutine  definition  statement  of  PTKSUC  is 
SUBROUTINE  PTKSUC  (KSUQ 

where 

KSUC  is  the  integer  variable  (1  <  KSX  <  NSX) 
defined  above. 

If  IPRINT  <  0  no  calls  are  made  to  the  output  subroutines. 

A  user  not  interested  in  the  use  of  any  one  of  the  output  subrou¬ 
tines  must  provide  the  corresponding  dummy  subroutine  (with  RETURN  as 
the  only  executable  statement)  in  order  to  avoid  unresolved  references 
problems. 
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FOPT  is  the  current  best  minimum  value  of  f  found  “from  algorithm 
start  (^Qp-p)  (FOPT  is  updated  whenever  a  function  value  f(x) 
is  computed) . 

FTFMIN,  FTFMAX  are  respectively  the  minimum  and  the  maximum  value 
of  f(x)  among  the  end-of-trial  values  obtained  at  the  final 
points  of  the  last  trajectories  of  the  current  trial 

fTFMIN^ ’ 

FTFOPT  is  current  minimum  value  of  FTFMIN  among  the  trials  which 
did  not  stop  due  to  the  stopping  condition  related  to  NPMAX 
(stopping  indicator  ISTOP  *  0,  see  below) .  FTFOPT  is  used  by 
SIGMA,  to  compute  the  input  parameter  KSUC  for  the  subroutines 
PTKSUC,  see  below, 

ISTOP  is  the  indicator  of  the  stopping  condition  of  the  trial,  as 
follows: 

ISTOP  *  0  The  maximun  nunber  NPMAX  of  observation  periods 
has  been  reached. 

ISTOP  ^  0  all  the  final  values  of  f(x)  of  the  last  obser¬ 
vation  period  (except  for  the  just  discarded  tra¬ 
jectory)  are  close  enough  to  their  common  minimum 
value  FPFMIN,  with  respect  to  an  absolute  or  rela¬ 
tive  difference  criterion,  ([1],  sect.  3.2.13), 
to  be  considered  nunerically  equal. 

If  ISTOP  f  0  the  absolute  value  and  the  sign  of  ISTOP  have  the 
following  meaning: 

The  absolute  value  indicates  which  of  the  difference  criteria 
was  satisfied 
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taken  place,  if  NSUC  (input  parameter  to  SIGMA)  had  been  given  a  (lower) 
value,  equal  to  the  current  KSUC. 


The  subroutine  PTKSUC  is  called  only  if  IPRINT  >  0  and  KSX  <  NSUC. 
The  subroutine  definition  statement  of  PTSEG  is 

SUBROUTINE  PTSEG  (N,  XPFMIN,  FPFMIN,  FPFNftX,  KP,  NFEV) 


where 

N  is  the  dimension  of  the  problem 

FPFMIN  and  FPFMAX  are  respectively  the  minimum  and  the  maximum  value  of 
f  (x)  among  the  values  obtained  at  the  final  points  of  the  trajec¬ 
tory  segments  of  the  current  observation  period  (excluding  the  dis¬ 
carded  trajectory) . 

XPFMIN  is  the  N-vector  containing  the  coordinates  of  the  final  point 
(or  possibly  one  of  the  points)  where  the  function  value  FPFMIN 
was  obtained. 

KP  is  the  total  number  of  elapsed  observation  periods  in  the  current 
trial. 

NFEV  is  the  total  number  of  function  evaluations  performed  from  algorithm 
start. 

The  subroutine  definition  statement  of  PTRIAL  is 


SUBROUTINE  PTRIAL  (N,  XOPT,  FOPT,  FTFMIN,  FTFMAX,  FTFOPT, 

ISTOP,  ISTOPT,  NFEV,  KP,  IPRINT) 

where 

N  is  the  dimension  of  the  problem 

XOPT  is  an  N-vector  containing  the  coordinates  of  the  point  (or  possibly 
one  of  the  points)  where  the  current  best  minimum  FOPT  was  obtained 


alleviate  the  efficiency  problems  connected  to  the  use  of  the  explicit 
Euler  step  on  ill-conditioned  functions. 

It  is  also  recommended  to  avoid  whenever  possible  to  provide  func¬ 
tions  such  that  the  "typical"  values  of  the  function  and  the  coordinates 
(rough  average  values  in  the  region  of  interest)  differ  from  unity  by  too 
many  orders  of  magnitude.  Such  a  care  is  generally  advisable  due  to  some 
numerical  values  adopted  in  the  FORTRAN  implementation,  for  example  to 
avoid  overflow,  but  may  be  absolutely  necessary  when  using  the  driver  sub¬ 
routine  SICMA1,  due  to  the  adopted  general  purpose  default  values  for  some 
input  data,  for  example  the  stopping  tolerances. 

3.5.2.  The  Output  Subroutines. 

Apart  from  the  output  parameters  in  the  call  statement  for  SICMA,  the 
package  is  designed  to  be  able  to  perform  external  output  also  by  means  of 
the  calls  to  three  output  subroutines  which  must  be  supplied  by  the  user: 
PTSEG,  PTRIAL,  and  PTKSUC.  The  calls  are  activated  according  to  the  value 
of  the  control  parameter  IPRINT  (sect.  3.2). 

The  subroutine  PTSEG  is  called  (if  IPRINT  >  0)  at  the  end  of  every 
observation  period. 

The  subroutine  PTRIAL  is  called  (if  IPRINT  >  0)  at  the  end  of  every 

trial . 

The  subroutine  PTKSUC  is  called  only  at  the  end  of  every  sucessful 
trial  such  that  an  increment  occurred  in  the  value  KSUC  of  the  maximum  num¬ 
ber  of  trials  which  had  a  uniform  stop  all  at  the  same  (current  or  past) 
value  of  ^Qpji  a  call  to  PTKSUC  therefore  provides  the  user  with  the  oper¬ 
ationally  interesting  information  that  a  final  success  claim  would  have 


3.S.  User-supplied  Subprograms. 

The  user  must  provide  the  function  FUNCT  which  must  compute  the 
value  of  f(x)  (sect.  1),  and  the  three  output  subroutines  PTSEG,  PTRIAL, 
PTKSUC.  The  above  subprograms  are  described  below:  all  non-integer 
arguments  are  double  precision  (integer  arguments  are  indicated  by  means 
of  the  FORTRAN  implicit  type  definition). 

3.5.1.  The  function  RJNCT 

FUNCT  nust  return  as  its  value  the  value  at  x  of  the  function 
f  to  be  minimized. 

The  function  definition  statement  is 

DOUBLE  PRECISION  FUNCTION  FUNCT  (N,X) 

where 

N  is  the  (input)  dimension  of  the  problem 

X  is  the  (input)  N-vector  containing  the  coordinates  of  the 
point  x  where  the  value  of  f  is  to  be  computed. 

Note  that  the  function  f(x)  should  comply  with  the  growth  conditions  (2.3), 
(2.4)  of  [1],  otherwise  the  fmction  must  be  suitably  modified;  this  may 
be  performed  by  simply  adding  a  penalization  function,  which  nust  be  zero 
on  the  region  of  interest.  We  note  that  this  device  can  be  used  also  to 
suitably  restrict  the  search  region  (for  example  in  the  case  of  periodic 
functions) . 

It  should  be  also  noted  that  —  although  some  form  of  automatic 
rescaling  is  provided  by  the  algorithm  —  it  is  certainly  advisable  to 
avoid  whenever  possible  to  provide  unnecessarily  ill-conditioned  functions 
(for  example,  due  to  careless  choice  of  physical  units),  in  order  to 


3.4  The  Driver  Subroutine  SIGMA1 


In  order  to  both  give  an  example  of  how  to  use  SIGMA,  and  to 
save  the  average  user  the  effort  of  deciding  the  numerical  values  for 
all  the  input  parameters  of  SIGMA,  a  driver  subroutine  SIGMA1  is  included 
in  the  package.  SIGMA1  sinply  calls  SIGMA  after  assigning  default  values 
to  a  number  of  input  parameters.  The  subroutine  definition  statement  is: 

SUBROUTINE  SIGMA1  (N,  X0,  NSUC,  IPRINT,  XMIN,  FMIN,  NFEV,  IOUT) 
where  the  parameters  have  the  same  meaning  as  in  SIGMA. 

All  the  other  input  parameters  of  SIGMA  are  assigned  default 
values  within  SIGMA1  as  follows: 


EPS  *  1 
DX  *  10"9 
IRAND  =  0 
NT RAJ  *  0 
ISEGBR  -  0 
KPBR0  »  0 
INKPBR  »  0 
NPMIN  «  10 
NPMAX0  *  100 
INPMAX  -  50 

NTRIAL  *  max(50,5*NSUC) 

TOLREL  -  10 “3 
TOLABS  «  10' 6 

KP.ASCA  -  10  (if  N<5);  -  300  (if  N  >  5) 

IXHP  -  1 

;  XRMIN(I)  *  -104  (I  *  1, . . . ,N) 

|  XRMAX(I)  «  104  (I  *  1 , . . .  ,N) 

i 

I 


NPMIN  (say  5  <  NPMIN  <  100) 

NPMAX0  (say  0  <  NPNIAX0  <  150) 

NTRIAL  (say  NTRIAL  >  50  and  NTRIAL  >  5-NSUC) 

INPMAX  (say  30  <  INPMAX  <  100) 

The  following  parameters  have  a  marked  effect  on  package  performance 
and  computational  effort: 

INHP,  NTRAJ,  ISEGBR,  INKPBR,  NSUC. 

The  magnitude  of  the  effect  roughly  decreases  from  left  to  right. 

In  order  to  avoid  untolerable  growth  of  the  computation  effort  or 
an  unacceptable  degradation  of  the  performance,  the  user  is  advised  to 
modify  (if  needed)  the  above  parameters  starting  from  XSUC,  based  on 
information  from  the  output  subroutines  (PTRIAL  and  PTKSUC) .  Note  that 
NSUC  is  the  only  "free"  control  parameter  of  the  driver  SI  (MAI  (Sect.  .4). 

The  value  of  KPASCA  should  be  based  on  possible  analytical  or 
experimental  evidence  on  the  ill-scaling  of  the  function  f(x) .  Choose  a 
small  value  (say  10)  for  a  badly  scaled  function,  a  large  value  (say  300) 
for  a  very  well  scaled  function.  The  N-dimensional  interval  (XRMIN,  XKMAX) 
should  be  as  large  as  possible,  consistently  with  the  purpose  of  avoiding 
computation  failures  (e.g.  overflow).  Finally  we  note  that  due  to  the 
joint  operation  of  the  stopping  conditions  for  the  trial  (see  [1] ,  sect. 
3.2.8),  in  order  to  use  only  one  of  the  conditions  it  is  sufficient  to 
put  to  zero  the  threshold  tolerance  (T0LREL  or  T0LABS)  of  the  other  con¬ 
dition.  Suggested  default  values  of  most  input  parameters  are  provided 
in  the  driver  subroutine  SI(MA1  (sect.  3.4). 


XMIN  is  an  output  N- vector  containing  the  coordinates  of  the  point 
(or  possibly  one  of  the  points)  where  the  final  value  FMIN  of 
fQpp  was  found. 

BUN  is  the  final  value  of  the  best  current  minimum  function  value 
fOPT* 

NFEV  is  the  (output)  total  nuirber  of  function  evaluations  (includ¬ 
ing  those  used  for  the  computation  of  derivatives,  and  for 
the  rejected  time- integration  steps). 

IOUT  is  the  (output)  indicator  of  the  stopping  conditions  as  follows: 
If  IOUT  *  -99  a  fatal  error  was  detected  when  performing  some  preliminary 
checking  of  the  input  data,  and  the  algorithm  was  not  even  started;  otherwise 
the  algorithm  was  started,  and  the  value  of  IOUT  is  the  final  value  of  the 
of  the  parameter  ISTOPT  (an  output  indicator  of  the  output  subroutine 
PTRIAL,  described  in  sect.  3.5.2.). 

Success  is  claimed  by  the  algorithm  it  IOUT  >0,  i.e.  if  at  least 
one  of  the  trials  stopped  with  a  positive  value  of  the  trial  stopping  indi¬ 
cator  ISTOP  (described  in  sect.  3.5.2)  and  no  lower  value  for  fQ^p  was 
found  in  the  following  trials . 

3.3.  Some  Guidelines  for  the  Choice  of  the  Input  Parameters. 

Proper  operation  of  the  package  should  be  almost  independent  of 
IRAND,  KPBR0  (and  X0) .  The  performance  of  the  package  should  not  be  too 
sensitive  to  H,  EPS,  DX,  since  these  are  initial  values  of  variables 
which  are  adaptively  controlled  by  the  program. 

The  following  parameters  are  expected  to  have  little  effect  on  the 
performance,  as  long  as  they  belong  to  wide  "insensitivity”  bounds: 


XRMIN,  XRMAX  are  input  N-vectors  defining  an  admissible  region  for  the 
x -values,  within  which  the  function  values  can  be  safely  com¬ 
puted  (see  [1],  sect.  3.2.11,  where  XRMIN (I),  XRMAX (I)  are 
caned  Rf,  RfV 

KPASCA  is  the  (input)  minimum  number  of  observation  periods,  before 
the  scaling  procedures  are  activated  (Kpasca) • 

IRAND  is  a  control  (input)  index  for  the  initialization  of  the 
random  nunber  generator: 

if  IRAND  >  0  the  generator  is  initialized  before  starting 
the  trial  Kt  with  seed  IRAND  +  1L-1; 

if  IRAND  <  0  the  generator  is  initialized  (with  seed  0)  only 
at  the  first  call  of  SICMA. 

INHP  is  used  to  control  the  number  NHP  ("duration”)  of  time  integra¬ 
tion  steps  for  observation  period  as  follows: 
if  INHP  ■  1,  NHP  ■  1  +  [log2(Kp)],  ("short"  duration) 
if  INHP  »  2,  NHP  ■  [v^  ]  ("medium"  duration) 
if  INHP  -  3,  NHP  -  Kp  ("long"  duration) , 
where  [x]  is  the  greatest  integer  not  greater  than  x. 

IPRINT  is  an  input  control  index  used  to  control  the  amount  of  printed 
output  by  controlling  the  calls  to  the  user-supplied  output  sub¬ 
routines  PTSEG  (end-of-segment  output),  PTRIAL  (end-of -trial  out¬ 
put)  ,  and  PTKSUC  (end-of-trial  output  related  to  the  count  of 
successful  trials),  described  in  sect.  3.5.2. 
if  IPRINT  <  0  no  call  to  the  print  subroutines 
if  IPRINT  -  0  call  only  PTRIAL  and  PTKSUC 
if  IPRINT  >  0  all  the  print  subroutines  are  called. 
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to  a  default  value  (NTRAJ  *  7),  and  if  the  input  value  is 
outside  the  interval  (3,20)  NTRAJ  is  set  to  the  nearest 
extreme  value) . 

ISEGBR,  KPBR0,  INKPBR  are  the  parameters  1^,  K^,  ^  which  determine 
which  one  of  the  simultaneous  trajectories  is  to  be  branched 
(see  [ 1] ,  sect .  3.2.6).  (Note  however  that  if  one  of  the  in¬ 
put  values  is  zero,  the  corresponding  variable  is  set  to  a 
default  value:  ISEGBR  *  (l+NTRAJ)/2,  (FORTRAN  integer  divi¬ 
sion)  ,  INKPBR  =  10,  KPBR0  *  3;  if  the  input  value  for  ISEGBR 
is  outside  the  interval  (1, NTRAJ),  ISE3GR  is  set  to  the  nearest 
extreme  value;  and  if  KPBR0  has  a  value  not  inside  the  interval 
(1, INKPBR),  it  is  assigned  the  same  value  modulo  INKPBR) . 

NFMIN  is  the  (input)  minimum  duration  of  a  trial,  i.e.  the  minimum 

number  of  observation  periods  before  checking  the  trial  stopping 
condition. 

NFMAX0  is  the  (input)  initial  value  (i.e.  for  the  first  trial)  for 
the  maximum  duration  of  a  trial,  i.e.  for  the  maximum  accept¬ 
able  number  NPMAX  of  observation  periods  in  a  trial  (N^^) . 

INPMAX  is  the  (input)  increment  for  NPMAX,  when  NPMAX  is  varied  from 
one  trial  to  the  following  one. 

NSUC  is  the  (input)  number  of  successful  trials  (with  the  same  final 
value  fgpp  see  sect.  1)  after  which  the  computation  is 
stopped  (Nsuc) . 

NTRIAL  is  the  (input)  maximum  allowed  number  of  trials,  after  which 
the  computation  is  stopped  (N^j^) . 

TOLREL  and  TOLABS  are  the  (input)  relative  and  absolute  tolerances  for 
stopping  a  single  trial  (x^,  x^,  see  [1],  sect.  3.2.8). 
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N,  X0,  H,  EPS,  DX, 

NTRAJ,  ISEGlX,  KPBR0,  INKPBR, 

NPMIN,  NPMAX0 ,  INPMAX, 

NSUC,  NTRIAL,  TOLREL,  TOLABS,  XRMIN,  XRMAX 
KPASCA,  IRAND,  INHP,  IPRINT 

SIGMA  returns  to  the  calling  program  the  output  parameters 
NTRAJ,  ISEGBR,  KPBR0,  INKPBR, 

XMIN,  BUN,  NFEV,  IOUT 

The  call  parameters  are  described  in  the  next  section. 

We  note  that  the  SIGMA  package  gives  the  user  the  possibility  of 
obtaining  —  during  algorithm  evolution  —  the  values  of  a  nunber  of 
other  parameters  by  means  of  the  output  subroutine  (to  be  supplied  by 
the  user)  which  are  described  in  sect.  3.5.2.  The  parameters  are 
KP,  NF,  XOPT,  FOPT 

XPBUN,  FPBIIN,  FPFMAX,  FTFMIN,  FTB1AX,  FTFOPT 
ISTOP,  ISTOPT,  KSUC 
and  are  described  in  sect.  3.5.2. 

3.2.  Description  of  the  parameters  of  the  call  statement  for  SIC^IA. 

N  is  the  problem  dimension  (nunber  of  coordinates  of  a  point  x) 

X0  is  an  N-vector  containing  the  initial  values  of  the  x-variables 
H  is  the  initial  value  of  the  time  integration  steplength. 

EPS  is  the  initial  value  of  the  noise  coefficient 

DX  initial  value  of  the  magnitude  of  the  discretization  increment 
(&x)  for  computing  the  finite -differences  derivatives. 

NTRAJ  is  the  nunber  of  simultaneous  trajectory  segments 

(Note  however  that  if  the  input  value  is  zero,  NTRAJ  is  set 
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3.  USAGE 


In  order  to  use  the  package  the  user  must  provide: 

a)  a  driver  program  which  calls  the  principal  subroutine  SIGMA, 

b)  a  set  of  four  auxiliary  subprograms  (to  compute  the  function 
f(x)  and  to  output  the  results) . 

The  CALL  statement  for  SIGMA  is  described  in  sect.  3 .1,  the  para¬ 
meters  of  the  CALL  statement,  are  described  in  sect.  3.2.  Some  guidelines 
for  the  choice  of  the  values  of  the  input  parameters  are  given  in  sect. 
3.3.  A  sample  driver  subroutine  (SIGMA1)  which  calls  SIGMA  is  described 
in  sect.  3.4:  such  a  subroutine  assigns  default  values  to  a  number  of 
input  parameters  to  SIGMA:  it  has  therefore  a  considerably  lower  number 
of  input  parameters,  and  can  be  used  as  an  easy-to-use  driver  for  the 
average  user.  The  user-supplied  subprograms  are  described  in  sect.  3.5. 

3.1.  Call  to  SIGMA 

The  call  statement  is 

CALL  SIGMA  (N,  X0,  H,  EPS,  DX, 

NTRAJ,  ISEGBR,  KPBR0,  INKPBR , 

NPMIN,  NPMAX0,  INPMAX, 

h'SUC,  NTRIAL,  TOLREL,  TOLABS,  XRMIN,  XRMAX, 

KPASCA,  I  RAND,  INHP,  I PRINT, 

XMIN,  FMIN,  NFEV,  IOUT) 

The  program  calling  SIGMA  must  set  the  input  call  parameters 
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